
////////////////////////////////////////////////////////////
// Routine to read a binary spincolor in the lime format
///////////////////////////////////////////////////////////


void read_binary_spincolor(complexspincolor *out, char *file)
{

  char *expected_lime_record;
  expected_lime_record = "scidac-binary-data";

  int nreals_per_spincolor=24;
  int nreals_tot=nreals_per_spincolor*L*L*L*T;

  int nbytes_per_spincolor_float=nreals_per_spincolor*sizeof(float);
  int nbytes_per_spincolor_double=nreals_per_spincolor*sizeof(double);
  int nbytes_per_spincolor_read;
  n_uint64_t read_bytes;

  LimeReader *reader;
  FILE *fp;
  int status;
  int read;

 //open the file
  fp = fopen(file, "r");
  if(fp == (FILE *)NULL) {
    fprintf(stderr,"Unable to open file %s for reading\n", file);
  }
  reader = limeCreateReader(fp);
  if( reader == (LimeReader *)NULL ) {
    fprintf(stderr, "Unable to open LimeReader\n");
  }

  read = 0;

  char *header=NULL;

  while( (status = limeReaderNextRecord(reader)) != LIME_EOF && read==0 ){

    if( status != LIME_SUCCESS ) {
      fprintf(stderr, "limeReaderNextRecord returned status = %d\n",
              status);
      exit(EXIT_FAILURE);
    }
    header=limeReaderType(reader);
    printf("found record: %s\n",header);
   if(strcmp(expected_lime_record,header)==0)
       {
          uint64_t nbytes=limeReaderBytes(reader);
          nbytes_per_spincolor_read=nbytes/(L*L*L*T);
          read_bytes=nbytes_per_spincolor_read/nreals_per_spincolor;
          if(nbytes_per_spincolor_read>nbytes_per_spincolor_double)
            {
              fprintf(stderr," ERROR The file contain %Ld bytes per site and it is supposed to contain not more than: %d or %d\n",(long long int)nbytes,nbytes_per_spincolor_read,nbytes_per_spincolor_double);
              break;
            }

           limeReaderReadData((char*)((double*)out),&nbytes,reader);
           read=1;

      }
}
  limeDestroyReader(reader);
  fclose(fp);


  if(nbytes_per_spincolor_read==nbytes_per_spincolor_float) //cast to double changing endianess if needed
    if(big_endian) floats_to_doubles_changing_endianess((double*)out,(float*)out,nreals_tot);
    else floats_to_doubles_same_endianess((double*)out,(float*)out,nreals_tot);
  else //swap the endianess if needed
    if(big_endian) doubles_to_doubles_changing_endianess((double*)out,(double*)out,nreals_tot);

}

void read_binary_colorcolorspinspin(colorcolorspinspin *S, char *filename_input)
{
  //loop over external spin-color
  char scfilename_input[1024];
  for (int is_ext=0; is_ext<4; is_ext++)  for (int ic_ext=0; ic_ext<3; ic_ext++)
  {
    int sc=3*is_ext+ic_ext;
    sprintf(scfilename_input,"%s.%02d",filename_input,sc);
    // Read spincolor in position space
    printf("Reading... %s\n",scfilename_input);
    complexspincolor *spinore=(complexspincolor*)malloc(sizeof(complexspincolor)*(T*L*L*L));
    read_binary_spincolor(spinore,scfilename_input);
    for (int t=0; t<T; t++)  for (int x=0; x<L; x++)   for (int y=0; y<L; y++)     for (int z=0; z<L; z++)
    {
       int sour=x+L*(y+L*(z+L*t)); // The spin color in position space is written in the input files with this convention (Francesco's convention)
       int dest=z+L*(y+L*(x+L*t)); // We will use this in the following to be compatible with tmLQCD
       for (int is=0; is<4; is++)   for (int ic=0; ic<3; ic++)
                S[dest][ic][ic_ext][is][is_ext]=spinore[sour][is][ic][0]+I*spinore[sour][is][ic][1];
    } 
   free(spinore);

  }
}
